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^\ . We review localization with non-Hermitian time evolution as applied to simple models of popu- 

lation biology with spatially varying growth profiles and convection. Convection leads to a constant 
[ imaginary vector potential in the Schrodinger-like operator which appears in linearized growth mod- 

els. We illustrate the basic ideas by reviewing how convection affects the evolution of a population 
influenced by a simple square well growth profile. Results from discrete lattice growth models in 
both one and two dimensions are presented. A set of similarity transformations which lead to exact 
results for the spectrum and winding numbers of eigenfunctions for random growth rates in one 
dimension is described in detail. We discuss the influence of boundary conditions, and argue that 
periodic boundary conditions lead to results which are in fact typical of a broad class of growth 
<—j ■ problems with convection. 
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03 . I. INTRODUCTION 

£ 

' Bacterial growth in a petri dish, the basic experiment of microbiology, is a familiar but interesting phenomenon. 
Depending on the nutrient and agar concentration, a variety of intriguing growth patterns have been observed |^-^|- 
O Some regimes can be modeled by diffusion limited aggregation, others by Eden models and still others exhibit ring 
O | structures or a two-dimensional modulation in the bacterial density. At high nutrient concentration and low agar 
density, there is a large regime of simple growth of a circular patch (after point inoculation), described by a Fisher 
' equation ||, and studied experimentally in Ref. g]. 

. Of course, most bacteria do not live in petri dishes, but rather in inhomogeneous environments characterized by, 
^sD ' e -g-i spatially varying growth rates and/or diffusion constants. Often, as in the soil after a rain storm (or in a sewage 
treatment plant), bacterial diffusion and growth are accompanied by convective drift in an aqueous medium through 
the disorder. By creating artificially modulated growth environments in petri dishes, one can begin to study how 
bacteria (and other species populations) grow in circumstances more typical of the real world. More generally, the 
challenges posed by combining inhomogeneous biological processes with various types of fluid flows seem likely to 

' attract considerable interest in the future. The easiest problem to study in the context of bacteria is to determine 
how fixed spatial inhomogeneities and convective flow affect the simple regime of Fisher equation growth mentioned 

' above. 

A derealization transition in inhomogeneous biological systems has recently been proposed, focusing on a single 
species continuous growth model, in which the population disperses via diffusion and convection the Fisher 
equation || for the population number density c(x, i), generalized to account for convection and an inhomogeneous 
growth rate, reads Q 

9c(x, t) /dt = DV 2 c(x, t) - v ■ Vc(x, t) 



+C/(x)c(x,i)-6c 2 (x,<), (1.1) 



where D is the diffusion constant of the system, v is the spatially homogeneous convection ( "wind" ) velocity, and b is 
& \ a phenomenological parameter responsible for the limiting of the concentration c(x, t) to some maximum saturation 
value. The growth rate U(x) is a random function which describes a spatially random nutrient concentration, or, for 
photosynthetic bacteria, an inhomogeneous illumination pattern ||. If £/(x) is constant over the entire sample, then 
the convection term — v • Vc(x, t) can be eliminated by a coordinate transformation and has no effect on the growth of 
the bacteria || . Only the introduction of a spatial dependence for the growth rate U (x) makes the convection term 
interesting. 

In this brief review, we summarize recent theoretical results concerning Eq. ( pTl| ). In Section [| we review results for 



a "square well" shape for the growth profile U(x) |10[. When linearized about the state of zero population, Eq. (1.1) 
is related to Schroedinger's equation with a square well potential and analytic results are possible. As we shall see, 
large enough convective velocities delocalize bound states in the equivalent quantum mechanical problem, and this 
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derealization is accompanied by escape of real growth eigenvalues into the complex plane. This import ant s ignature of 
derealization is preserved for random growth profiles. A lattice approximation to the continuum Eq. (1.1), useful for 
treating both square well and random growth profiles, is described in Section III. Exact results and various similarity 
transformations for the lattice model in one dimension Jll|l are summarized in Section IV. The one dimensional lattice 
model is used to explore the effect of boundary conditions in Section In the appendix we give analytical expressions 
for the spectrum of a two dimensional homogeneous system for two different lattices and the continuum model. 



II. RESULTS FOR SQUARE WELL GROWTH PROFILES 



In the following we review results for equation Eq. (1.1) with a simple "square well" growth profile U(x.), imposing 
a positive growth rate a on an illuminated patch ("oasis"), and a negative growth rate — ea outside ("desert") flCf] : 



C/(x) = 



for Ixl < f , 



— ea , for Ixl > -5- , 



(2.1) 



where W is the diameter of the oasis. Experimentally this situation can be realized using a very simple setup, 
which illustrates the basic ideas of localization and derealization and leads to interesting further questions. A one 
dimensional example is shown in figure [jj where a solution with photosynthetic bacteria in a thin circular pipe, or 
annular petri dish, is illuminated by a fixed uniform light source through a mask, leading to a "square well" intensity 
distribution. The mask is moved at a small, constant velocity around the sample to simulate convective flow. (Moving 
the mask is equivalent to introducing convective flow in the system, up to a change of reference frame ||.) The bacteria 
are assumed to divide in the brightly illuminated area ("oasis") at a certain rate, but division ceases or proceeds at a 
greatly reduced rate in the darker region ("desert") outside. As a result, the growth rate in this continuum population 
dynamics model is positive in the oasis and small (positive or negative) in the surrounding desert region. Using this 
simple nonlinear growth model, one can study the total number of bacteria expected to survive in the steady state, 
the shape of their distribution in space and other quantities, as a function of the "convection velocity" of the light 
source llQ]. 




FIG. 1. Experimental setup: a solution with photosynthetic bacteria in a circular pipe or a thin annular track in a petri 
dish is illuminated only in a small area, while the rest of the sample is either kept dark or illuminated with reduced intensity. 
The light source (shaded rectangle) is moved slowly around the sample to model convective flow. The bacteria are assumed 
to divide in the illuminated area ("oasis") at a certain growth rate a > 0, and die (or grow modestly) in the remaining area 
("desert") with growth rate — ea. 
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Much can already be learned by linearizing equation (1.1) around c(x, t) = 0: 



dc(x,t)/dt = £c(x,t), (2.2) 

with the linearized growth operator 

£ = £>V 2 -v-V + ?7(x). (2.3) 

The convection velocity v acts formally like an imaginary vector potential in this Schr odinger-like equation jl2| . For v 
nonzero, C is non-Hermitian, but it can still be diagonalized by a complete set of right and left eigenvectors {^(x)} 
and {0^(x)}, with eigenvalues r n , and orthogonality condition 

J d d x4,i(x)4>X(x)=5 m , n , (2.4) 

(d is the dimension of the substrate, we focus here on d = 1 or d = 2). The time evolution of c(x, t) is then given by 

c(x, t)=J2 c «^(x) exp(r„i) , (2.5) 

n 

where the initial conditions and left cigenfunctions determine the coefficients {c n }, 

d d a;^(x)c(x,t = 0) . (2.6) 



(See Refs. |8[] and 10 for a discussion of how the results from the linearized problem can be used to treat the full 
nonlinear system.) 



Fig. H shows the complex eigenvalue spectrum associated with Eq. (2.3) for a one-dimensional system with the 



potential (|]l]) and periodic boundary conditions, at four different values of the convection velocity v. In Ref. JTcj ] 
details of the derivation of these results are given. At zero velocity, C is Hermitian and all eigenvalues r„ are real. 
There are bound states (discrete spectrum) and extended or delocalized states (continuous spectrum). Because C 
resembles the negative of a quantum mechanical Hamiltonian, bound states have the highest eigenvalues (growth 
rates), in contrast to the usual situation in quantum mechanics. At finite velocities, all delocalized states, except 



the uppermost, acquire a complex eigenvalue. According to Eq. (2.5), states with positive real part of the eigenvalue 
(ReT n > 0) grow exponentially with time, states with negative real part (ReT n < 0) decrease exponentially with time. 
In a large one dimensional system the "mobility edge" |13j , which we define to be the eigenvalue of the fastest growing 
delocalized state, (i.e. the rightmost eigenvalue in the complex parabolas of figure |J), is located for one dimensional 
systems of size L at the overall average growth rate 

r* = (U)=[ dxU{x)/L ~ -ea, (2.7) 
Jo 

In Fig. ||, the eigenvalues of the localized states compose the discrete, real spectrum to the right of the mobility 
edge. With increasing velocity these localized eigenvalues move uniformly to the left by an amount v 2 /AD, and 
successively enter the continuous delocalized spectrum through the mobility edge, which remains fixed. The parabola 
also broadens in the vertical direction: the imaginary parts of the eigenvalues of the delocalized states grow by an 
amount proportional to v. A given localized right cigcnfunction <j>^ undergoes a "derealization transition" when the 
velocity reaches a corresponding critical derealization velocity v — i>*, at which its eigenvalue T n has been shifted so 
far to the left that it just touches T*. At higher velocities it joins the parabola of eigenvalues describing a continuum 
of delocalized states. The ease with which such a derealization transition can be observed experimentally depends 
on whether there are growing delocalized eigenstates in the system, i.e. whether the mobility edge has a positive real 
value or not. 
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Spectrum for different drift velocities 
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FIG. 2. Complex nonhermitian eigenvalue spectra (normalized by the difference of the growth rates inside and outside the 
oasis Uo = a + ea — 1) at velocities above and below the extinction transition. The spectra are extracted from numerical 
simulations of a lattice approximation to the continuum described in Section |ll]|, for a one dimensional system of 1000 sites. 
In units such that the lattice constant £q = 1, we have oasis width W = 20 sites, diffusion constant D — 0.3Uo£q, with growth 
rate — ea — — 0.5t/o in the desert and a = +0.5Uo in the oasis, so that the average growth rate is — 0.48t/o (which is equal to the 
mobility edge T* up to finite size effects). The velocity parameters (in units of the Fisher wave velocity in the oasis vp = 2\faD 
||) are v/vp ~ (circles), v/vp ~ 0.39 (crosses), v/vp ~ 0.93 (triangles), and v/vp ~ 1.39 (squares). The point spectra 
are slightly offset in the vertical direction so as to be able to distinguish the eigenvalues of the localized states for different 
velocities. As the convection velocity increases, the mobility edge remains fixed, and the parabola of the delocalized eigenvalues 
opens up as v/vp is increased. The real, localized eigenvalues move to the left for higher velocities. When v/vp > 1.39 all 
states are delocalized. 

We distinguish three cases, (U) < 0, (U) > 0, and (U) ~ 0: 

(1) In a large "deadly" desert (J.U) ~ — ea < 0) all delocalized states die out, because the mobility edge lies to 
the left of the origin, as in Fig. The growth rate of each localized eigenstate 4>n then becomes negative at a 
corresponding "extinction" velocity v nc which is smaller than the corresponding derealization velocity u*. Thus, as 
convection is increased, the population dies out before it can delocalize. Total extinction occurs when the eigenvalue 
of the localized "ground state" (fastest growing eigenfunction of C) passes through the origin, i. e. when v = vq c (with 
VOc > v nc for all n > 0. The total bacterial population of the steady state in the corresponding nonlinear problem 
goes to zero linearly with v — vq c Jl0| . 

(2) If the average growth rate (U) is positive (i.e for a small enough desert or a small positive growth rate in 
an infinite desert), the mobility edge lies to the right of the origin and the delocalization transition can indeed be 
observed at v — Vq where the "ground state" becomes delocalized. One expects to see universal behavior near this 
delocalization transition, since there is a diverging correlation length in the system, which renders microscopic details 
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irrelevant for certain quantities. For example, the correlation length £ of the spatial density distribution scales as 
£ ~ l/(v — VqY with a universal exponent v = 1 fl(i[| . 

(3) A special (universal) behavior is expected for the spatially average growth rate (U) = 0. In this case the 
derealization and extinction velocities coincide. Fig. [| summarizes the different scenarios in a sketch of the phase 
diagram for large systems with fixed well depth Uq = a + ea, obtained by tuning the drift velocity, and the average 
growth rate (U) ~ — ea. Also shown in Fig. || is a horizontal transition line at (U) — separating a small velocity 
region (e < 0) where localized modes dominate the steady state bacterial population, from one (e > 0) containing a 
mixture of localized and extended states. It is of course also possible to drive a population extinct at zero velocity 
simply by lowering the average growth rate. This special transition at (U) = — U c is indicated at the bottom of figure 
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FIG. 3. Schematic phase diagram in one dimension for infinite system size, as a function of average growth rate (U) and 
convection velocity v for fixed "well depth" Uo = a + ea. For a deep well (Uq S> D/W 2 ) the extinction transition out of the 
localized phase occurs when v = for (U) = —U c , where U c — Uo- The diagram shows that if the growth rate is negative 
outside and inside the oasis (i.e. (U) < —U c ), then the only possible state is extinction at any velocity. If there is positive 
growth inside the oasis, but negative outside, a localized population can survive in the oasis, but only for small enough wind 
velocities v. Extended states are present for a small positive growth rate in the desert ((U) > 0). In this case localized and 
delocalized states coexist for small velocities ("mixed phase"), while at large velocities all eigenstates are extended. The ground 
state becomes delocalized at the critical velocity v c which marks the phase boundary between the mixed and the extended 
phase. 



III. LATTICE MODEL OF GROWTH AND CONVECTION 

As discussed in Refs. M and [juJ, analytic expressions for eigenf unctions and eigenvalues obtained from an analysis 
of a linearized growth problem can be used to study the nonlinear steady state provided only a few unst able growing 
modes are present. In more general situations, however, a numerical analysis of the full nonlinear Eq. ([hi]) may be 
necessary. Numerical work also provides insights into situations where the growth profile is random, as opposed to 
the simple square well considered in Section y|. We review here a lattice approximation to the continuum nonlinear 
model which is very useful for obtaining numerical results. It was used in Ref. JlCfl to obtain growth spectra like those 
shown in figure ^ to confirm and extend analytic results obtained in one dimension. Here, we apply the lattice model 
to some illustrative two dimensional growth problems with convection. 

A discrete lattice approximation (discrete in space but continuous in time) to Eq. (1.1) in d dimensions reads j|] 
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^1 = w J2[eZ- e »c x+e „(t) + e-^c^Ct) - 2cosh(g ■ e u )c x (t)] 

+t/(x)c x (i) - bcl (t) , (3.1) 

where c x (i) is the species population at the sites {x} of a hypercubic lattice with lattice constant £q, and the {e^} 
are unit lattice vectors. Furthermore, w ~ D/£q, where D is the diffusion constant of the corresponding continuum 
model, and g ~ — v£q/(2D), where v is the conv ective flow rate of the continuum model. U(x) and b have the same 



interpretation as in the continuum model (1.1). The subtraction in the first term insures that c x (i) is conserved 
c x(0 == 0) if U(x) =6 = 0. When linearized about c x = 0, Eq. (3.1) may be written 



^eM = ££(x,x')c x ,(i) (3.2) 
where the discrete Liouville operator C is the matrix 

d 

1 = W ^5>-s-e,| x + e ,)( x | + e g ' e "|x)<x + e!/ |] 

X V—\ 

+ £L4 + t/(x)]|x)(x| (3.3) 

X 

with A = — 2wJ2t=i cos h{g. ■ e„). Typical spectra with a square well growth profile for a one dimensional 1000-site 
model were shown in Fig. ]^ for several values of g oc v |Io|] . Parameter values and energies studied were chosen to 
give a good approximation to the continuum. 

For two dimensional systems, a triangular lattice (see Eq. (|3.5| )) is actually a better approximation to the continuum 
than a square lattice, because it minimizes artifacts caused by preferred lattice directions. Using the unit lattice vectors 
e i, e 2, e 3 of Fig. H to define the Laplacian, one finds Jl4| via a Taylor's series expansion 

2 3 



=i 



= V 2 c x (0 + ^(V 4 c x (t)) + 0$) . (3.4) 

Using the lattice vectors ei, and e2 of a square lattice gives a similar result, but the fourth order terms in that case 
are £q(8 a /dx A + <9 4 /<9?/ 4 )c x (t)/12 which cannot be expressed in terms of the Laplacian of c x (i) [Q. 
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FIG. 4. Hexagonal oasis (filled circles) surrounded by a desert (empty circles), in 2 dimensions on a triangular lattice, with 
periodic boundary conditions. The inset in the lower left corner of the figure displays the three basis vectors e„ of the lattice. 
The arrow in the upper right corner of the figure indicates the 30° direction of the convection velocity v relative to the horizontal 
axis. (This direction renders a better approximation to the continuum than for example a direction along one of the three 
lattice vectors e„.) 

Fig. H shows the eigenvalue spectra at finite convection velocity v for an oasis in a desert for two different lattices, 
for t he s quare lattice, with v along the x axis(a), and along the 45°-axis (b), and for the triangular lattice (see 
Eq. ( |3.5| )) with v along the x axis (c) and along the 30°-axis (d) as in Fig |[ States localized near the oasis again 
appear as a string of real eigenvalues above a pattern of complex eigenvalues representing delocalized eigenfunctions. 
Although the ellipsoidal pattern of extended state eigenvalues is an artifact of the lattice, the upper edge can be a 
good approximation to the continuum if the direction is chosen appropriately. Directions coinciding with "lattice 
planes" lead to spurious features in the delocalized part of the spectrum. Spectral artifacts are particularly noticeable 
for the square lattice, when v is along one of the lattice directions (a). Note also, that for the square lattice the 
ellipse of complex eigenvalues is mirror symmetric with respect to the two coordinate axes through the center of the 
ellipse. For neither the triangular lattice, nor the continuum model, there is such a mirror symmetry with respect to 
the imaginary axis in the complex spectrum. An easy check is given by an analytic computation of the spectrum for 
a homogeneous growth rate in the system. The results are cited in the appendix. 
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Im{r}/U versus Re{r}/U for square and triangular lattices 
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FIG. 5. Complex nonhermitian eigenvalue spectra for a two dimensional oasis in a desert, using different lattices. The model 
parameters (with £o = 1) are D = w = 1, v/vf = 0.37, a = 29.5w, and ea — 0.5w. (a) and (b) show results for a square lattice 
with 400 sites and a square oasis with 36 sites, with periodic boundary conditions. In (a) the velocity v is directed along the 
a;-axis, which is one of the two preferred lattice directions. In (b) the velocity is directed along the diagonal (45°-axis). (c) and 
(d) show the results for a triangular lattice with 400 sites, a hexagonal oasis with 37 sites, with periodic boundary conditions, 
as shown in Fig ^. In (c) the velocity is directed along the x axis, which coincides with the direction of one of the lattice vectors 
e, y . In (d) the velocity is directed along the 30°-axis, as indicated in Fig (The eigenvalues are normalized by the difference 
of the growth rates inside and outside the oasis Uo = a + ea — 30w.) 

We now use a tri angular lattice to study the derealization transition in two dimensions. For a triangular lattice 
the sum in Eq. (3.3) over v is replaced by a sum over three unit lattice vectors illustrated in Fig || One obtains 



|x + e„)(x| + e*- e "|x)(x + e„| 



+ 5> + £/(x)]|x)(x| 



(3.5) 



with A = — ^ Ylu=i cosh(g-e u ). Fig. ^ shows lattice simulation results for six spectra at different convection velocities 
for the linearized two dimensional system with a hexagonal oasis, and periodic boundary conditions. The convection 
velocity is oriented along the 30° axis relative to the lattice, as shown in figure 0. 
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Spectra for an oasis in 2 dimensions 

( lm{rj/U vs Re(r}/U for 6 different convection velocities ) 
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FIG. 6. Complex nonhermitian eigenvalue spectra at six different drift velocities for a two dimensional system with a hexag- 
onal oasis (with growth rate a = 29. 5w', where w' = 1 and w' = 2w/3) in a desert (with negative growth rate — ea = — 0.5w') 
with periodic boundary conditions. The growth rates of the spectra are divided by the difference of the growth rates inside 
and outside the oasis Uo = a + ea = 30w'. The drift velocities are directed along the 30° axis of the triangular unit cell. Their 
respective values (divided by the Fisher wave velocity inside the oasis vf = 2y/aD) are indicated next to each figure. The 
spectra are extracted from numerical simulations of the triangular lattice model with 400 sites, with a hexagonal oasis with 37 
sites, and diffusion constant D — 1.5w'£q, with £o = 1. (This figure actually only shows the part of the spectrum which gives a 
good approximation to the continuum problem.) 

At zero velocity all eigenvalues are real (figure [|(a)). Just as in one dimension, with increasing convection velocity, 
the real growth rates of the localized states decrease while the mobility edge (i.e. the rightmost eigenvalue of 
the parabolic envelope of the delocalized states with complex eigenvalues) remains unchanged (b). Unlike in the 
one dimensional case, there are delocalized states with similar real parts of their eigenvalues, but a spread in the 
imaginary parts. (The reason is that in two and higher dimensions there are delocalized states that can carry a 
quantum mechanical current in directions other than the convection velocity JI^JTsfl .) Also, as the convection velocity 
is increased, some of the localized states become delocalized (signaled by eigenvalues which escape into the complex 







plane jlj]) even before joining the parabola of delocalized states (c). The parabola formed by the eigenvalues of 
the delocalized states widens as the imaginary parts of the eigenvalues grow with the convection velocity, just as we 
saw for one dimensional systems. For even higher convection velocity the eigenvalues of the localized states enter 
the parabola of the delocalized states. The mobility edge is ill defined at these intermediate velocities (d). As the 
convection velocity v approaches the Fisher wave velocity vf = the upper edge of the band settles again into 

its fixed value (approximately Re{T} — for the given parameter values), where it remains for all higher velocities 
(f)- 

Dynamical properties of the full nonlinear system can be studied by discretizing Eq. (3.2) also in time and using 
Runge Kutta methods to solve the equation numerically [|l6 17 l||,[l8|]. Although we hope to present a more thorough 
study at a later date, Fig. ^ shows preliminary results in the form of snapshots of a bacterial colony drifting and 
expanding in a 2 dimensional system on a triangular mesh, with convection and quenched random growth rates. The 
initial condition for the system shown was a delta function like spatial distribution concentrated in the upper right 
hand corner. Theoretical predictions in the limit of high convection velocity for the long time scaling behavior of a 
growing and drifting droplet of bacteria in the linearized regime were given in [0 . 
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FIG. 7. Contour plots for snapshots of the distribution of bacteria in a two dimensional landscape with (quenched) random 
growth rates, obtained from the full nonlinear problem for convection in the direction relative to a triangular lattice shown 
in Fig. Q with velocity parameter g = 1, w' = 2w/3 = 1, b = 1, and io = 1. The random growth rates U(x) are uniformly 
distributed with values ranging from —0.3w' to 0.7w' . At time t=0 the bacterial population was zero everywhere, except at one 
lattice point where it was set to 1. The snapshots show the expanding and drifting droplet after 250, 500, 750, and 1000 time 
steps of size 5t — 0.1/w'. The contour lines marking population density c = 0.1 are the innermost contour for the snapshot 
after 250 time steps, the second to outermost contour for the snapshots after 500 and 750 time steps and the outermost contour 
for the snapshot after 1000 time steps. The simulated system had 500 x 500 sites. 



IV. EXACT RESULTS FOR LATTICE MODELS IN ONE DIMENSION 



Consider the lattice models discussed in the previous section with a random growth rate U ( x) chosen from an 
identical proability distribution at every lattice site x. If the drift parameter g oc v in Eq. (|3.3| ) is set to zero, the 
growth operator C becomes formally identical to a Hermitian Hamiltonian describing a quantum mechanical particle 
which hops between lattice sites in a disordered medium. A number of exact results are known for this problem in one 
dimension; we review here recent exact results for the non-Hermitian generalization embodied in Eq. ( |3.3| ) Jl9|-pl|, 
with an emphasis on the winding numbers and complex currents associated with extended states [jTT| . 

A numerical study of the eigenvalue spectrum of the lattice model in one dimension |l^] reveals the following: when 
g = 0, all eigcnfunctions ip n {x) are localized with a bounded point-like spectrum {r„} of eigenvalues lying on the real 
axis, similar to the subset of localized states which fall to the right of the mobility edge for the square well spectrum 
shown in Figj2[ Localization of all states by a random potential is the expected result for quantum particles in a 
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disordered medium in one dimension. The situation changes dramatically, however, when g 0. For g greater than 
a critical value g = g c , a bubble of complex eigenvalues appears in the center of the band, similar to the response 
to convection for arbitrary g of the continuous spectrum par t of the spectrum for a square well shown in Fig.|[ The 
eigenf unctions associated with this bubble are extended |12| . Eigenvalues outside this bubble remain real, and are 
related by a simple rigid band shift to their values for g = 0. The right and left eigenfunctions <f>^' L (x) for states 
outside the bubble remain localized, and are related to the eigenfunctions ip n (xj) of the Hcrmitian problem for g = 
by 



,±9 x i 



(4.1) 



We assume here a one dimensional lattice with unit spacing £q = 1. 

We shall focus on the eigenfunctions and complex currents associated with the band of extended states in one 
dimension for lattice growth models with disorder in one dimension [pd| . Unlike delocalized states in Hermitian disor- 
dered systems (where the eigenfunctions can always be chosen to be real), we show that these complex eigenfunctions 
are characterized by a conserved winding number n, even when the disorder is strong. Such topological quantum 
numbers can be used to label the eigenvalue spectrum T n (g), where g is the asymmetry parameter. A study of the 



eigenvalue trajectories as a function of g then leads to complex currents defined by J n = — 
response of the n-th eigenvalue to changes in the drift velocity. 

The matrix representation of the one dimensional non-Hermitian Liouville operator (| 
localized at positions {xj , j = N} reads, 
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(4.2) 



where we assume periodic boundary conditions, and the constant background of the growth rate A has been omitted. 

The matr ix displayed in Eq. (4.2) (tridiagonal, with nonzero corner matrix elements) is actually a generalization 
of Eq. (3J3): As before, Uj is a zero mean random growth rate chosen independently for each site. However, the 



hopping rate between sites i and i + 1, «;,, is now a random function of position, which allows for fluctuations in the 
local diffusion constant. We have also allowed for fluctuations in the local convection velocity gi. This latter variation 
can be eliminated via similarity transformation which preserves the eigenvalues of C, and produces only a smooth, 
nonsingular change in the eigenfunctions. It is helpful to carry out this transformation in two stages. The first isolates 
all effects of convection on the bond connecting the iVth and 1st site: Upon defining 
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we have 
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Except for the corner matrix elements, £ describes Hermitian population growth without convection. A second 
similarity transformation, 
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where g = X^7=i then redistributes the convection uniformly on all bonds, 



The composite similarity matrix 
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(4.6) 



Sij = (S 1 S 2 )i ] = %exp[-y^ff fc 



fc=i 



A' 



(4.7) 



fc=i 



thus leads to a population growth problem whose eigenvalues reflect a mean convection parameter g = ^X^Lj.fl'.r 
For the g/s chosen independently at each site, the sample to sample fluctuations should fall off like a/N where a is the 
variance of the distribution. As N grows one may thus replace the fluctuating quantity gi by a disorder independent 
average value g, as we shall do in the rest of this paper. This result implies universality in the response to random 
convection in one dimensional growth models - nonuniform convection velocities may b e m apped into a uniform 
average velocity via a similarity transformation, provided one applies the transformation (4.7) to the eigenfunctions 
as well. 

For a ID ring with random parameters {Uj\, {wi} and g = 0, all the eigenfunctions of ( |4.4| ) are real and localized, 
and its eigenvalues are real and discrete p3,||. We assume for simplicity that all Wi > and that the chain is large 
but finite, such that although the spectrum is discrete, the length of the chain is much larger than the maximal 
localization length of an eigenmode. 

In a typical one-dimensional disordered system with g = 0, the localization length £ (defined via the exponential 
fall off of the eigenfunction) is larger at the center o f th e band and smaller at the tails. The criterion for derealization 
of the asymmetric system via convection (see Eq. (4.1)) is k£q < <?, where k = l/£ |l^] . As a result, pairs of complex 
energies representing delocalized states first appear as a "bubble" at the center of the band, which then spreads into 
the band tails as the convection parameter g is increased (see Fig. ||) . To study the complex currents associated with 
delocalized states we follow Refs. fl9|-|2~l|l and exploit the relation between the complex spectrum of the asymmetric 
problem with g 7^ and the real eigenvalues of a "background" localized problem with g — 0, with the same realization 
of disorder. 
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FIG. 8. Eigenvalue trajectories in the complex plane. Points on the real axis are eigenvalues for the Hermitian problem 
for g = . The complex eigenvalue t rajec tories T n (g) and F* (g) are orthogonal to the level curves (dashed lines) denned by 
Eq.(4.1C). Inset: Angles entering Eq.( [l.9[ ) for a given eigenvalue T(g). 



The condition for a complex number T = Tr + iTi to be an eigenvalue of the matrix C(g) displayed in Eq. (4.6) 
is Det[TI — C(g)] = 0, where / is the identity matrix. Upon expanding in cofactors, we find that T is an eigenvalue 
provided pIS 



N 



Det[FI-L(g = 0)} = l[(r-r o i ) 
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N 

2[cosh( 5 iV)-l]n( ! |) 



(4. 



where the {T®} are the (real) eigenvalue s of the background matrix C(g — 0). To extract winding numbers, we first 
observe that the rig ht h and side of Eq. (4^) is real and positive. As a result, the phases of the complex numbers in 
the left product of 



for each complex T should sum up to (see inset to Fig. 

r«-r? 



pir, 



(4.9) 



with p = 2n, where n is an inte ger, and the function cot~ 1 (x) varies from ir to as x goes between — oo and oo. 

As Ti — ► 0, each term in Eq.(4.9) gives ir for every eigenvalue T° to the right of Tr, and zero for each eigenvalue 
to the left (see Fig.|J). To satisfy (4.9), the eigenvalue must leave the real axis and enter the complex plane at the 
gap between the 2n-th and the 2n+l-th eigenvalues of the g = "background" system. We call n the index of the 
trajectory Tji(g) + iTj(g) in t he c omplex plane. 

We see immediately from (4.9) that the rightmost eigenvalue (with n = 0) must remain real, consistent with 
Perron- Frobenius theorem |23| . The corresponding nodeless eigenfunction corresponds to the ground state of —C 
For N even, particle-hole symmetry implies that both the rightmost and l eftm ost eigenvalues are always real. More 
generally, for a fixed value of n, the set of all [Tfj(g), Tj(g)] satisfying Eq. (4.S) defines a curve in the complex plane, 
as illustrated in Fig.^. Henceforth, we assume N even for simplicity. 

A more complete description of the eigenvalue trajectory results from taking the logarithm of the modulus of Eq. 
(4.8). In the limit Ng >> 1 one finds a second constraint on T(g), namely, 



ln(- 



) + ^EM|r-r° 



(4.10) 



where |r - = \J (Tr — T 1 /) 2 + T'j and w = 2 Qli=i IT F- 

This constraint is described graphically in Fig. |8J, which shows schematically the level curves defined by Eq. (|4~To| ) 
near the band center for three values of g. These are lines of constant potential for an equivalent 2d electrostatic 
problem with charges at the positions of the localized point spectrum for g — 0. When g is small, the constraint is 
solved by eigenvalue pairs on the real axis in the gaps between neighboring r°'s indexed by p = 2n. As g increases, 
successive pairs of eigenvalues eventually merge at a saddle point in the potential contours and detach from the real 
axis at right angle. For a given real energy r^, Thoulcss has defined an energy dependent inverse localization length 
K(Tft) for the associated g = Hermitian problem p2j, 



«(r fl ) = -in(|) + l^in(|r fl -r°|) 



(4.11) 



Upon comparing with Eq. (4.10), we see that Ti becomes nonzero whenever 
ment point of the eigenvalue pair. 



g\ > k(T' r ) where (r' fl ,0) is the detach- 



As g increases above k(T'), Eqs. (4.9) and (4.10) define a unique pair of complex eigenvalue trajectories T n (g) and 
* (g) for every value of n. Upon passing to t he lim it N — > oo, the spectrum for g = closes up, and is described 

_ Jft and ( pl| ) may then be recombined into a single complex equation, namely 



by a density of states po(X). Eqs. (p 



dX po(X) ln[r n (p) - A] = In 



(?") 


■ f 2n \ 


+ «7T 







(4.12) 



where 



ln[r(<7) - A] = ln[\T R + iT T - A|] + icor^Ta - A)/rj] 



(4.13) 



For simplicity we set the average growth rate to zero. Then in the limit N » 1 and for a density of states function 
symmetric around the special detachment point with T' R = 0, there is a purely imaginary trajectory of the form 
T(g) = iTi(g) with n = N/2. For fixed g, Eq. Q4.12[ ) thus leads to an implicit formula for the "height" Tf ax of the 
bubble of complex eigenvalues in the center of the band (see Fig. m 1 namely, 



dX p(X) ln[A 2 



Mf). 



(4.14) 
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This integral vanishes, as expected p(J, in the "one way" limit, g — > oo with ^e' 9 ' = 1. F or ot her detachment points, 
the eigenvalue trajectories curve to the left or right as required by the constraint of Eq. ( 4.1 0| ) (see Fig. 



The analysis above suggests that the imaginary parts of all {T n (g)} (except the two which remain on the real 
axis), diverge as \g\ — * oo. It is then expected that all eigenfunctions 4> n (j) are approximately plane wav es, <j> n (j) ~ 
exp(ik n rj /£q), with free particle eigenvalue energies T n (g) =wcos{k n + ig) |,[l§. For large {T^g)], Eq. Q leads 
to 

w 

T n (g)^ 1 exp(\g\+2i7rn/N). (4.15) 

Comparison with the free particle spectrum at large \g\ shows immediately that the index n of the eigenvalue trajectory 
and the wave vector are related, k n = 2im/N . These wave eigenfunctions spiral around the origin in the complex 
plane as one moves along the Id lattice of the tight binding model sites, leading to a well defined winding number n. 
The winding number associated with the eigenvalue r* (g) in the lower half plane is then —n. 

As g decreases, it can be shown that the associated delocalized wave function <j> n {j;g) must remain nonzero for 
all j, so long as the eigenvalue retains an imaginary part p4j . Thus, the winding number must be preserved as g 
decreases, i.e., the winding number is a topological invariant along an eigenvalue trajectory. The projection of such a 
delocalized eigenfunction is illustrated in the inset to Fig. 0. 

As N — > oo, we can replace the winding index n by the continuous variable k n — jrn. Eq. ( |4. 1 2| ) then shows quite 
generally that the complex spectrum is a function only of the combination g + ik, i.e., T n (g) = T(g + ik). It then 
follows from the Cauchy-Riemann relations that all eigenvalue trajectories are at right angles to the lines of constant 

dF(g) -drfg) 




FIG. 9. Spectrum and winding numbers for TV = 1000 with Uj uniformly distributed in the interval [— 1, 1], and g — 0.4. 
Inset: Projection onto the complex plane of an extended eigenfunction near the center of the bubble with N = 200, g = 1, and 
winding number n = 100. 



An explicit formula for the complex current results, moreover, from differentiating Eq. ( 4.12 ) with respect to g 



; dT n (g) 
dg 



A) (A) 



= -iG^iTnig)} 



(4.16) 



Evidently, the current associated with a particular complex eigenvalue T n (g) is determined by the Green's function 
of the 5 = problem, evaluated at the complex energy T n (g). 

Fig. ^| shows a typical eigenvalue spectrum for C in one dimension superimposed on the winding numbers associated 
with the extended eigenfunctions. The inset illustrates how the winding number is defined for a typical extended 
state. The winding numbers are exactly ±N/2 at the band center, and their magnitudes decrease monotonically as 
one moves toward the upper edge of the band, i.e., toward the lowest energies of the corresponding Hamiltonian. 
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The winding numbers remain finite up to the mobility edge separating complex and real eigenvalues, but become 
undefined in the band tails, where all wave functions are real. At the transition, an overall constant phase factor can 
be adjusted so that the imaginary part of the eigenvector vanishes while its real part remains finite. In this case, the 
winding index classification can be replaced by counting the nodes of the real, localized wave function ]22] . 



V. PSEUDOSPECTRA AND SENSITIVITY TO BOUNDARY CONDITIONS 



In this review, we have explored the spectral properties of the linearized growth operator C in order to learn 
about the time evolution of the growing populations using spectral decomposition, as in Eq. ( |2.5| ). In the recent 
years, however, it has been recognized t hat such a straightforward analysis of non-Hermitian operators and matrices 
(such as the Liouville operator in Eq. (|2.3|)) can in some circumstances lead to inaccurate predictions [p5|-p7|. In 
this section we discuss convection-diffusion with a random growth rate in this context. We argue that, while some 
non-Hermitian problems arc indeed abnormally sensitive to small perturbations, our results (provided the boundary 
conditions admit a net population flux) are quite robust. We consider a linear, (possibly non-Hermitian), operator 
which admits matrix representation in a Hilbert space. When convenient, we interchange freely between the matrix 
and the operator descriptions. 

Consider first the Hermitian operator, 7i, which has a complete set of orthogonal eigenvectors. If the spectrum of 
H is non-degenerate, any eigenvalue T n corresponds uniquely to an eigenstate. When a two-fold degeneracy occurs, 
there are two independent eigenvectors corresponding to the same eigenvalue, and one may orthogonalize these states 
using the Graham- Schmidt scheme |23]. If Ti is an evolution operator, its response to external perturbation with 



"frequency" z is determined by the resolvent {i.e. the Green's function) (zl — 7i) _1 . The resolvent norm |28| diverges 
if z is a spectral point of TL. For small deviations, i.e., z = T n + Sz, for some r„ one may use the diagonalized 
representation of the resolvent to see that its norm is proportional to (Sz) . 

For a non-Hermitian operator, or, in particular, a real asymmetric object, such as C(g), the situation is different. 
If the spectrum of this non-normal operator is non-degenerate, it admits two complete, biorthogonal states of left and 
right eigenfunctions For any eigenvalue T n (g) the left eigenstate < c/>£| satisfies < 02|£(g) = T n (g) < </>2| while 
the right state satisfies £(g)\(f)^ i > = r„(g)|<^, >. The biorthogonality condition means that 

< CI 05 >= <W- (5-1) 

Unlike the Hermitian case, the left and right eigenstates are not simple complex conjugates. Hence, the normalization 
condition (|]l]) implies no restriction on the usual Hilbert norm of the left and right eigenfunction considered separately, 
and the only normalized quantity is the "overlap integral" between them. In the following we explore the significance 
of this observation. 

Non- normality also implies that, in the generic case, as two eigenstates coincide (when the spectrum is "degenerate") 
the corresponding eigenvectors also coincide. The subspace spanned by these vectors is one dimensional, and an N x N 
matrix is defective, i.e., it fails to span the N dimensional vector space associated with it. An illustration of the "almost 
degenerate" situation for the left and right eigenvectors in a 2 x 2 space is repr esented in Fig. |l^. Note that < _Z>2 1 is 



orthogonal to \Ri >, and < L x | is orthogonal to \R 2 > as required by Eq. (5.1) 
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FIG. 10. Two pairs of almost degenerate, left and right eigenvectors of a non- normal matrix. The condition (5.1) implies 
that the Ri and L2 are orthogonal, as well as R2 and L\. Close to degeneracy Li and L2 almost coincide, hence Ri is "almost 
orthogonal" to L\. As illustrated in Fig. [y], a typical value of the poperly normalized eigenfunctions must then be very large 
in order to satisfy the normalization condition < L n \R n >= 1. 



The notion of pseudospectra has been introduced p5| in order to quantify the problems which may arise when 
spectral analysis is applied to non-normal operators. While the spectrum T(£) of the operator C is defined as the set 
of points (in the complex plane) which are eigenvalues of C, the e-pseudospectrum T e (C) is defined as a set of points 
which imitate a true eigenvalue to within a small parameter e > 0. Formally, the definition is: 



r e (£) = {complex numbers z such that \\(zl — C) || > e } 



(5.2) 



) 1 = 00). Eq. ( pT2| ) implies that the response of the system to a small time 
z mimics a resonance up to times of order t ~ 1/e. If the pseudospectral 



(note that if z is a spectral point, (||z/ — C\ 
dependent perturbation with "frequency" 
region is large for small e, the system is highly unstable against weak external perturbations even far away from the 
resonance. 



Another definition (mathematically equivalent to (5.2)) is, 



r e (£) = {complex numbers z which belong to the spectrum of C 
for all 5C such that \ \SC\ \ < e}. 



sr., 



(5.3) 



Eq. (5.3) defines the pseudospectrum via the sensitivity of C to small perturbations in the growth operator. Again, 
broad pseudospectral regions for small e imply high sensitivity to small perturbations. 
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region 

FIG. 11. Sketch of the right and left eigenfunctions of a non-normal operator. Since the normalization is given by the overlap 
(shaded) region, the spectrum of such a system is very unstable against a non-diagonal perturbation which connects the regions 
X and Y. 



The reason for such an extreme sensitivity in the context of the random growth operator considered here, can be 
clarified if one looks at the (localized) left and right states shown in Fig. [0]. For simplicity, we show the "ground 
state" (most rapidly growing) left and right eigenfunctions, but similar considerations apply to excited states as well. 
For nonzero convection velocities the peaks in the eigenfunction will be shifted away from the overlap region according 



to Eq. (4.1). If the overlap integral denned in Eq. fc.1) is small, the amplitude after normalization of both 02 011 the 
left side of Fig. [ll] and of (fy at the right side of Fig. ll| will be very large. Perturbation theory will then yield huge 
corrections to the eigenvalues if 5L has an off diagonal matrix element < cf>L\5C\(j)R > connecting regions X and Y. 
The pseudospectral instability thus becomes stronger and stronger as the overlap integral < 4>™\4> r R > becomes small 
compared with the norm of the individual eigenvectors. 

Fig. |l0| illustrates why this overlap integral also becomes smaller as eigenvalues approach degeneracy. According 
to the biorthogonality condition ([O]), if n ^ m, the right eigenvector </>^ should be orthogonal to all the left states 
0£ . Near the degeneracy point, when the two left eigenvectors almost coincide, the right eigenvectors become nearly 
orthogonal to their own left counterparts, and the overlap integral almost vanishes. Thus, nearly degenerate eigenvalues 
in a non-Hcrmitian operator can cause extreme sensitivity to small off-diagonal perturbations, corresponding to large 
pseudospectral regions for small e. The same mechanism may have caused the unstable spectra of large non-Hermitian 
matrices (with random entries connecting sites arbitrarily far apart) observed recently by Chalker and Mehlig pc| ]. 

With this background in mind, let us consider the pseudospectral properties of the convection-diffusion operator 
with random growth considered in this review. The pseudospectrum of the convection-diffusion operator 

C cd = 2 x + vd x (5.4) 

has been considered by Reddy and Trefethen ^?J, for "absorbing" boundary conditions. In our language this situation 
corresponds to a uniform growth rate with no population allowed outside of the [0, d) region. As for an infinite 
potential well in quantum mechanics, the eigenfunctions satisfy (f)(0) = 4>(d) = 0. For v = 0, C c d is Hcrmitian with 
a complete set of orthogonal eigenfunctions, indexed by an integer m, d > m ~ sin(rmrx / d) with the corresponding 
eigenvalues T m = —ir 2 m 2 /d 2 . As v is increased, the transformation ( |4.l| ) gives left and right eigenfunctions of the 
form <fim R ~ exp(±ua;) sin(rn7ra;/(i), leading to peaks near the boundaries, similar to those sketched in Fig. [n]. The 
spectrum undergoes a rigid shift of — v 2 /A for all m, since the above transformation always works for these boundary 
conditions. The result of this left-right asymmetry is, as discussed above, a broad pseudospectral region. In fact, as 
shown in Ref. [ p7| (for the case v = 1), all points in the complex plane contained inside the parabola 

Re z = -(Im z) 2 (5.5) 



are contained in T e (C c d) for small e, and the norm in Eq. ( |5.2| ) diverges exponentially as d — > oo. In the light of our 
discussion above - this extreme instability arises from small perturbations of the evolution operator which connect 
the two bumps at the edges, i.e., from small changes of the boundary conditions. 

We show that the results in this review are in fact unaffected by this interesting instability. First, we note that for 



periodic boundary conditions, the eigenstates of the non-Hermitian operator C c d (see Eq. (5.4)) are simply 4>^' L = 
exp(±ik n x) with k n = 2irn/d and eigenvalues T n — —k^ + ivk n , i.e., the eigenvalues lie on the parabolic limit of 



the pseudospectral region illustra ted in Eq. (5.5) for v = 1. Since the wave functions in this case do not have the 



asymmetric form implied by Eq. (4.1), we expect stability against small perturbations. 
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Consider the random matrix model of Sec. IV, now with a change of a single Wj to allow boundary c ond itions 



similar to those of Ref . |27| . Specifically, we insert a factor 77 in front of the corner matrix elements of Eq. 

wn — > i] w N ■ (5-6) 
By tuning 77 to zero, we break the chain and discard the periodic boundary conditions. Upon making the same 



transformation which led to Eq. (4.4), we see that all effects of convection vanish when 77 = 0, so the eigenvalue 



spectrum remains real for an arbitrary array of convection velocities gi , similar to the results for "absorbing" boundary 



conditions in the simple convectio n-di ffusion problem in Ref. |27|. If 77 is set to zero in Eq. (4.4), it is easy to show 



that a 'gauge' transformation like (4.1) relates the eigenfunctions to their g = values for all values of g, with a trivial 



shift of the real 3 = eigenvalue spectrum. 

The results change dramatically, however, for any 77 > 0, no matter how small. Indeed, it is clear from Eq. (4.4) 
that the eigenvalue spectrum for nonzero 77 is similar to that for 77 = 1, provided we calculate the spectrum with 
effective values of g ± , 

g ± = g±ln(r))/N, (5.7) 

where the + sign refers to the lower left corner matrix element and the — sign refers to the upper right. All effects of 
77 vanish in the continuum limit N — > 00, and the eigenvalues approach a universal spectrum like that in Fig. ^|, for 



all 77 ^ 0. It is straightforward to show, using the methods of Sec. fV, that the eigenfunctions are also unchanged in 
the large N limit. We conclude that the results for periodic boundary conditions are stable to a change in boundary 
conditions in one dimension, provided the chain is not actually broken by setting 77 = 0. 

Our conclusions about the stability of the spectrum and eigenfunctions for periodic boundary conditions also apply 
to the simple convection diffusion operator without randomness. The exceptional sensitivity to small perturbations of 
the spectrum for "absorbing" boundary conditions p7[ can be understood in the following way: Since the operators 
we consider here are real, all eigenvalues must be real, or occur in complex conjugate pairs. If, in addition, we impose 
absorbing boundary conditions, we exclude all states carrying a real current. This exclusion implies real eigenfunctions 
which insures that all imaginary parts of the eigenvalues vanish for the lattice model in the limit 77 = 0. Replacing 
this absorbing boundary condition by a weak link allows the eigenvalues to escape into the complex plane. The real 
spectrum studied in pTj , considered as the 77 — )• limit of a lattice model, involves the simultaneous coincidence of all 
eigenvalues in pairs (except those at the band edge). This condition of extreme defectivity (an N x N matrix with 
only ~ N/2 linearly independent eigenvectors !) is responsible for the exponential instabilities studied in p7J. 
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APPENDIX A: ANALYTICAL RESULTS FOR THE SPECTRA IN THE HOMOGENEOUS SYSTEM 

We give here the expressions for the spectra of a 2 dimensional system with a homogeneous growth rate, computed 
for a square lattice, a triangular lattice, and the continuum model. For a homogeneous growth rate the eigenstates 
can be assumed to be plane waves, with k-vectors in the first Brillouin zone. For simplicity we set the growth rate 
everywhere to zero. (A nonzero constant growth rate merely amounts to a constant shift of the eigenvalues T(k).) 
For the square lattice one obtains 

2 

T(k) = 2w [cosh(g ■ e„)(cos(k • e„) — 1) + isinh(g ■ e v ) sin(k ■ e v )\ . (Al) 

For the triangular lattice the eigenvalues are 

Aw 3 

r(k) = — y^[cosh(g • e„)(cos(k • e v ) - 1) + i sinh(g • e„) sin(k • e„)] , (A2) 
while for the continuum model the expression is simply 
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T(k) = -Dk 2 + ik • v , (A3) 
a form which agrees with the expansions of Eq. ( Al ) and Eq. (|A2|) for small wave vectors k. 
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